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We study instabilities of single-species fermionic atoms in the p-orbital bands in two-dimensional 
optical lattices at noninteger filling against interactions. Charge density wave and orbital density 
wave orders with stripe or checkerboard patterns are found for attractive and repulsive interactions, 
respectively. The superfiuid phase, usually expected of attractively interacting fermions, is strongly 
suppressed. We also use field theory to analyze the possible phase-transitions from orbital stripe 
order to liquid-crystal phases and obtain the phase diagram. The condition of nearly-perfect Fermi- 
surface nesting, which is key to the above results, is shown robustly independent of fermion fillings in 
such p-orbital systems, and the {2kF, ±2fcF) momentum of density wave oscillation is highly tunable. 
Such remarkable features show the promise of making those exotic orbital phases, which are of broad 
interest in condensed-matter physics, experimentally realizable with optical lattice gases. 
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I. INTRODUCTION 



II. SYSTEM AND MODEL 



Ultracold atoms in optical lattices are highly tun- 
able quantum systems, which are ideal to simulate con- 
ventional condensed-matter physics. Traditional stud- 
ies carried out on optical lattices only involve the low- 
est s band. In recent years, studies on higher orbital 
bands of optical lattices have shown many interesting re- 
sults For bosons, staggered orbital orders in square 
lattices and stripe orders in triangle lattices [1] have 
been proposed. Most recently, experimental progress 
has realized superfluidity on higher-orbital bands [7H9|, 
and experimental signatures distinguishing the staggered 
orbital order have been proposed in recent theoretical 
work 0, [13, [IH- For fermions, frustrated orbital or- 
ders in the strongly interacting Mott regime were found 
years ago. Recently, novel Fulde-Ferrell-Larkin- 
Ovchinnikov (FFLO) phases [11 [H on the p-orbital 
bands [l6l [TtI and multiband superconductivity [Tsj in 
optical lattices have been proposed, keeping this field fas- 
cinating. 

In the present work, we study the interacting spinless 
p-orbital fermionic atoms in two-dimensional (2D) square 
optical lattices with both attractive and repulsive inter- 
actions. We find that the quasi-one-dimensional feature 
of the Fermi surfaces of the double degenerate Px- and 
Pj,-orbital bands gives rise to the following interesting or- 
ders. For attractive interactions, it induces charge den- 
sity wave (CDW) order in a wide filling regime where the 
superfiuid order is greatly suppressed. For repulsive in- 
teractions, orbital density wave (ODW) order is induced. 
Both CDW and ODW show stripe or checkerboard pat- 
terns in space, depending on the filling. We further show 
that our system is a simple, clean, and highly tunable 
system to realize possible nematic and smectic liquid- 
crystal phases, which is a topic of great current interest 
in correlated condensed-matter physics [l9l - [2^ . 



Consider a system of spinless fermions filled up to de- 
generate Px- and Py- orbital bands in a 2D square lat- 
tice. Such a system can be realized by considering an 
anisotropic three-dimensional optical lattice with lattice 
potential Kp = ^u=x,y,z^vSvc?{kLr^), where fc^ is the 
wave vector of the laser beams and the lattice constant 
is a = 7r/fcL. By setting Vz ^ Vx = Vy, we realize dy- 
namically decoupled 2D square-lattice layers in the xy 
plane, each being a 2D system. The 2D system is then 
filled with spinless fermions such that the lowest s band 
is fully occupied and two degenerate Px- and Py- orbital 
bands are partially filled. In general, the band gap be- 
tween the s and p bands is much larger than the interac- 
tion, and the s-band fermions are dynamically inert. By 
expanding the fermionic field operators in the Wannier 
basis and using a tight-binding approximation, we obtain 
the p-band Fermi-Hubbard model 

ra/3 

-ti^na,r + g^nx^rny^r (1) 
ra r 

to describe the system with chemical potential ^. Equa- 
tion ([1]) only contains nearest-neighbor hopping and on- 
site interaction, since in typical ultracold-atom experi- 
ments the next-nearest-neighbor hopping and nearest- 
neighbor interaction are negligible. In Equation ((1}, 
Ca^r is the annihilation operator of Wannier state Pa at 
site r, and ria^r = C*! r^a,r is the number operator for 
the Pa-orbital state at site r. The subscripts a and /3 
run over x and y. The hopping term tajS is given by 
tap = ['^ll'Ja/S — t±{^ — ^ap)] , whcrc the parallel (trans- 
verse) hopping <|| {t±) means the hopping of pa-orbital 
fermions at site r to the nearest neighbor r -f e^, with 
P — a {(3 ^ a). Here, Gq is the lattice unit vector in 
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FIG. 1. (Color online) A schematic diagram illustrating how 
the {2kF ,2kF) momentum of density fluctuation satisfies the 
nesting Fermi-surface condition. Here the lattice constant a 
is set to be one. Red (dark gray) solid curve: Fermi surfaces 
of pi,-orbitaI band. Green (light gray) solid curve: Fermi sur- 
faces of py-orbital band. Blue dashed line: Fermi momenta 
of px- and Py- orbital bands. Solid arrow: the {2kF, ^kp) mo- 
mentum of density fluctuation simultaneously satisfying the 
nesting Fermi-surface condition for both p^- and Py- orbital 
bands. 



the a direction. The last term is the onsite interaction 
between p^- and Py-orbital fermions induced by p-wave 
scattering, with g as the coupling constant. Due to the 
p-orbital anisotropy, we expect i|| ^ In quantum 
chemistry, t|| {t±) is referred to as the a (tt) bond. Wc 
can use a harmonic approximation to estimate t|| and 
t±, where a standard tight-binding calculation gives the 
transverse hopping tj_ = e^^''^^-' Vx/2, and the parallel 
hopping i|| = I'T]'^ /2 — l\t±. The parameter 77 = ax a is 
typically a large number (^ 1) and therefore i|| 3> t±. 
Here, ai, = {V^/ Ep!)^/^kL, where Er = h?kj^/2m is the 
recoil energy for atoms of mass m. The onsite interac- 
tion is given hy g = gpalaz{22al -I- a^)/32(27r)^/^ in the 
pseudopotential approach with coupling constant gp [l^l . 
In the paper, the lattice constant, Boltzmann constant, 
and Planck constant are all set to be one. 



III. FERMI SURFACE INSTABILITIES 

It is well known that nesting Fermi surfaces are crucial 
to realize some spontaneously translational-symmetry- 
brcaking phases, e.g., CDW, spin density wave (SDW), 
and FFLO. For the lowest s band in a 2D square lattice, 
the nesting Fermi surfaces only occur at half filling, as- 
suming only nearest-neighbor hopping. In contrast, in 
our system, the nesting of quasi-onc-dimensional Px and 
Py Fermi surfaces, as shown in Figure [1] is independent 
of filling for a wide range of /i, as long as t± <C In 
Figure [TJ px and Py Fermi surfaces are perpendicular to 
each other, which greatly suppresses the Cooper insta- 
bility from particle-particle channel scattering. The rea- 
son is that in order to induce Cooper instability, all of 
the fcrmion pairs need to have almost the same center- 
of-mass momentum, which is impossible here with each 



particle-particle pair composed of one Px- and one py- 
orbital fermion, given only onsite interaction in Equa- 
tion dD). In contrast, eachp^, (py) particle- hole pair in the 
density channel is composed of one particle and one hole 
within the Px- (Py-) orbital band, which benefits from the 
nesting Fermi-surface condition. To simultaneously sat- 
isfy the nesting Fermi-surface condition for both px- and 
Py- orbital bands in the density channel, the momentum 
of density fluctuation should be 



Ql,2 ~ {2kF,±2kF), 



(2) 



as shown by the black arrow in Figure [1] where kp is the 
Fermi momentum for each band. 

To illustrate the above statement quantitatively, we 
study different instabilities in our system by random- 
phase approximation (RPA) as follows. The validity of 
RPA is discussed in Appendix [X] For the density chan- 
nel, we define the density operator p^.q = '^I fc+gC'a.fc 
in momentum-Matsubara frequency space for the p^- or- 
bital band where a ~ x,y. The 2-1-1 momenta k and 
q are defined as fc = (k, ia;„i) and q = (q, iWn), where 
= (2m -|- l)7rT and a;„ = 2n7rT are fermionic and 
bosonic Matsubara frequencies. The density-density cor- 
relation function in our system without interaction has 
the form 



T 



k+q 



(3) 



where {■■■)^ means thermal average without interaction. 
The spectrum of Pa;-orbital fermions is (^x,]s. — 2t|| cos kx — 
2t± cos ky — fi, and ^j,.k has a similar form. The Fermi 
distribution function np is given by np{S^k) = 1/(6^*=/"^ -I- 
1). At q = Q = {2kF,±2kF), t_L (perfect nesting) 
and iiJn = (static limit). Equation ^ reduces to 
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(4) 



in the continuous limit. In Equation D{Ep) is the 
density of states near the Fermi surface and ujd is some 
energy cutoff. 

When the interaction is turned on, the density-density 
correlation function IIq^ (g) ~ -^{pa,qPi3,-q) can be eval- 
uated by RPA, as shown in Figured! where (...) means 
thermal average with interaction. We then evaluate the 
following correlation functions: 



T 

n± = -^{p±{q)p±{-q)), 



(5) 



where p± [q) = px^q^Py^q are the total density and density 
difference between p^; and Py fermions, i.e., the CDW and 
ODW instability channels. Given that ^ 0, q = Q = 
{2kp, ±2kp), and icun ~ 0, Equation ^ reduces to 
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2X° 



(6) 
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FIG. 2. (Color online) Feynman diagrams of density-density 
correlation functions of the same orbital band lixx and differ- 
ent orbital bands Hxy, and the pair-pair correlation function 
Dpp . Red (dark gray) and green (light gray) lines are the 
propagators of free px and py fermions. 



Since 



X" ^ Dln(^j, any arbitrarily small attractive 
(repulsive) interaction g < (g > 0) can induce diver- 
gence of n+ (n_) at sufficiently low temperature. Such 
divergence indicates phase-transition to the correspond- 
ing symmetry-breaking phase. 

In experiments, a small but finite t_L is inevitable, 
which makes the Fermi-surface nesting not perfect. The 
density wave orders do not exist even at T = if the in- 
teraction is too small. However, for a nonperfect nesting, 
the same density wave orders are generally expected to 
exist if the interaction strength exceeds a certain critical 
value. For example, even if the Fermi surfaces between 
spin-up and -down fermions are Zeeman split in the pres- 
ence of a field, it is well known that the BCS superflu- 
idity or superconductivity persists up to a critical Zee- 
man splitting for a given interaction strength. The latter 
is known as the Chandrasekhar-Clogston limit [23l [2^. 
Increasing \g\ is experimentally feasible due to the high 
tunability of the interaction in optical lattices, e.g., by 
increasing the lattice potential. 

Figure [3] shows the phase-transition temperatures from 
CDW instability evaluated by RPA with small transverse 
hoppings, t± = 0,0.04, and 0.08. We set the system size 
to be N'^ — 300^. The parallel hopping is set to be 
i|| = 1 as the energy unit, and we choose the interaction 
strength 5 = — 2. It can be seen that although a small 
t± = 0.04 weakens the (stripe) density wave order, at a 
finite interaction g = —2, the (stripe) density wave order 
still occurs. However, a stronger t± = 0.08 destroys the 
(stripe) density wave orders over a certain range of the 
effective chemical potential (including Hartree term) fi' . 

The particle-particle (Cooper) channel can be studied 
in a similar way by evaluating the correlation function of 
the pair operator, Ag = J2k^x,-k+qCy,k- The pair-pair 
correlation function in our system without interaction 
reads 

n°p('z) = |i(AtA,)", 

L 1 ~ " F(Ca;,-k+q) " 'T-F(?y,k) 



(7) 
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FIG. 3. (Color online) phase-transition temperature for CDW 
instability at (j = —2 from RPA calculation. t± — (a) 0, 
(b) 0.04, (c) 0.08. Red (dark gray) line: instability towards 
the checkerboard density wave order. Green (light gray) line: 
instability towards stripe density wave order. The transi- 
tion temperature towards checkerboard density wave order 
with the effective chemical potential (including Hartree term) 
fi' near is much higher than that towards the stripe den- 
sity wave order with fi' away from 0, which indicates that 
the former is much stronger than the latter. This feature 
comes from the Umklapp process at half filling. Besides, the 
phase-transition temperature towards the checkerboard den- 
sity wave order does not show any noticeable change as t± in- 
creases from to 0.08, which suggests that the checkerboard 
density wave order is not affected by t±. 



without interaction, by choosing q = Q = {2kp,zL2kp), 
Equation ^ has logarithmic divergence as shown in 
Equation (jj]). In contrast, no logarithmic divergence is 
found in Equation ([7]) at any value of q. Therefore, the 
pair-pair correlation with interaction. 
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pp 



(8) 



evaluated by RPA as shown in Figure [H will not di- 
verge at any temperature, which is different from the 
density-density correlation functions in Equation ([6|). It 
means that the instability of the particle-particle chan- 
nel is greatly suppressed, and there is no phase-transition 
towards superconductivity. 



IV. MEAN FIELD THEORY AT T = 

The above consideration only shows that a phase- 
transition towards density wave order can happen in our 
system. In order to find the ground-state property, i.e., 
the order parameter, we apply a real-space mean-field 
analysis at T = for both g > and 5 < 0. The interac- 
tion part of the Hamiltonian given by Equation ([Ij can 
be decoupled in the density channel such that 



Recall that for a density-density correlation function 



where {na,r) = A^a.r is the self-consistent condition and 
(...) means the expectation value of the ground state at 
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FIG. 4. (Color online) The spatial total density pattern of 
CDW obtained from real space mean-field analysis (showing 
12^ out of 300^ sites) with p = -2 at (a) = and (b) = 
— 1. For ODW {g = 2), the density difference shows similar 
patterns to (a) and (b) for /x = and n = 1, respectively, (c) 
A schematic LDA phase diagram in the presence of a trap. 
|r| is the distance from the center of the trap. I, II, and 
III are the regions of stripe, checkerboard, and stripe density 
wave orders, respectively. Here x, y, and r represent the site 
numbers. 

T = 0. Eq. (m reduces to 

rafi ra 

+.9^(?^x,r^^y,r + ny^y-M^^^ - M^^y-My^r), 
r 

(10) 

which is in quadratic form and can be solved self- 
consistently. We set the parameters the same as be- 
fore in Sec. IIIH with = to simplify the calcula- 
tion. For attractive interaction g = —2, we find CDW 
order where the densities of Px and Py fermions are the 
same. When ^ = 0, the total density of the ground 
state exhibits a stripe pattern in real space as shown in 
Figure IDJa), and the energy per site is —2.0295, lower 
than that of the homogeneous-density state of —2.0282. 
When ^ = —1, the total density exhibits a checker- 
board pattern as shown in Figure IHJb), and the sys- 
tem is at half filling. In this case, the ground-state en- 
ergy per site is —0.7826, which is lower than that of 
the homogeneous-density state of —0.7731. Fourier se- 
ries Ma,r = o.o + J2'^=i i'^n cos(nq • r) + bn sin(nq ■ r)] are 
then used to fit Figure |H[a) and (b). It can be seen that 
the chemical potential is modified to fi' = fi — a^g by the 
background density ag (Hartree term), and the filling is 
determined by ^' instead of ^. We find q « (0.427r, 0.427r) 
by fitting FigurelU^a) . Higher-order harmonics (the n > 1 
Fourier components) are found to be nonvanishing in this 
case as expected in CDW, but are very weak compared 
to the first-order terms (a„ ^ ai, 6„ <C fei, n > 1). 
The checkerboard pattern in Figure SJ^b) has momentum 
q = (tt, tt). The q's in both cases agree with Equation ([2]) 
very well, with kp determined by the effective chemical 
potential n' . 

For repulsive interaction g = 2, we find ODW order, 
where the densities of Px- and Py-orbital fermions are no 
longer the same, and the difference between them (ODW 
order parameter) oscillates in space. When /i = and 
/i = 1, the density difference shows stripe and checker- 
board patterns similar to Figure m^a) and (b). 



There are some general properties of the density wave 
orders. (1) The checkerboard order at half filling is much 
stronger than the stripe order with other fillings due to 
the Umklapp process, which greatly enhances density 
wave order at half filling [2^. (2) In general, increas- 
ing transverse hopping t± weakens the nesting Fermi 
surface condition and tends to destroy the stripe order. 
However, the checkerboard order at half filling is not af- 
fected by the Fermi-surface curvature. It is because a 
(tt, tt) momentum of checkerboard order always satisfies 
the perfect nesting condition, independent of t± [25j . The 
above two features can be studied by calculating instabil- 
ities as shown in Sec. Illll or evaluating energy gain with 
mean-field analysis. (3) The Neel orbital order found in 
Ref. [m at half filling and strong-coupling limit g — > +oo, 
where Px and py Wannier orbitals alternate in space, can 
be understood as the extreme case of checkerboard ODW. 

In experiments, a shallow harmonic trap V{r) is 
present in addition to the optical lattice, and a spatial 
phase separation is expected due to the additional trap- 
ping potential. With local density approximation (LDA) 
/i(r) ^> /i(r) — V{r), a schematic phase diagram is shown 
in Fig. m^c). At the center region of the trap where the 
local chemical potential is the highest, the stripe order 
exists due to large filling (region I). As one moves to- 
wards the edge of the trap, the filling decreases, and 
when the effective local chemical potential fJ,'{r) « 0, 
the checkerboard order appears and this region is at half 
filling (region II). As one moves towards the edge fur- 
ther, the filling becomes low and the stripe order emerges 
again (region HI). Therefore, our theory predicts a spa- 
tial density profile of phase separation with stripe core 
— checkerboard shell — ?• stripe edge. 



V. LIQUID-CRYSTAL PHASES AT T / 

The above mean-field analysis shows the existence of 
density wave orders at T = 0. As one raises the tempera- 
ture, the thermal melting effect may drive the system to 
different liquid-crystal phases before it becomes normal 
Fermi liquid. In the following part, we will first present 
a field theory which incorporates the thermal melting ef- 
fect to study the liquid-crystal phases in square-lattice 
systems. Then we will make connections between this 
field theory and the specific microscopic model discussed 
before, i.e., how to determine the coefficients in the field 
theory from Equation ([T]). At last, we will comment 
on the advantages of our system to study liquid-crystal 
phases. 

For simplicity, we only consider liquid-crystal phases 
from stripe CDW with attractive interaction, where the 
densities of Px- and py-orbital fermions are the same. 
ODW with repulsive interactions are left for future study. 
The stripe CDW breaks the C4 rotational symmetry of 
the square lattice down to C2 (a Z2 phase-transition), 
and also breaks lattice translational symmetry, as shown 
in Figure |3Ja). As a result, two types of topological de- 
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fects can occur at finite temperature: the Z2 domain 
walls and the (edge) dislocations of stripes, where the 
latter may drive the system to smectic or nematic liquid- 
crystal phases. The smectic liquid-crystal phase breaks 
both translational symmetry and C4 rotational symme- 
try, which is essentially the same as the stripe order, 
while the nematic liquid-crystal phase only breaks C4 ro- 
tational symmetry and can be viewed as melted smectic 
stripes ^ [13 . 

The total density fluctuations associated with 
momenta Q1.2 can be parameterized as Sp = 
[(/iie**^^ "" -I- 02e*'^^ '" -I- C.C.] . For the incommensurate 
case, by adopting Gaussian approximation and keeping 
only up to quartic terms, the effective action reads 

+V\(f>i\^\<f>2\^ + St + (11) 

where St denotes topological defects of stripe disloca- 
tions, similar to the vortex term in the XY model. 
Eq. (jlip is invariant under C4 rotations of 7r/2, tt, and 
37r/2, which yields ((?!)i, 02) (02, 0i), (</>t , 02)7 and 
(02, 0i), respectively. Besides, Eq. pT|) also has two U(l) 
symmetries, i.e., 0i^2 01,26*'''^'^, where ipi and (p2 are 
arbitrary global phases. The coupling constants j, r, u, 
and V can be derived from the microscopic model given 
by Equation ([T]) as shown in Appendix [Bl and we find 
V = 4u, which strongly suppresses the coexistence of 0i 
and 02 . Without loss of generality, we assume 02 is sup- 
pressed, i.e., the saddle point is at |0i| ~ $ and |02| ~ 0, 
and write 0i = ^e'^'^ . Neglecting the (gapped) amplitude 
fluctuation of $, the low-energy theory is described by ip 
as S^ = i/ci2r.7$2(v^)2. 

The smectic and nematic order parameters can be 
defined as (0i — 02) and (|0ip — I02p), respectively. 
At T = 0, (01 - 02) - (01 ) - $e''^« ^ 0, with ipo 
as an arbitrary global phase, and the system is smec- 
tic. At arbitrary small temperature, the gaplcss U(l) 
mode of ip restores the translational symmetry, causing 
(0i) = (02) = 0. The system is algebraic smectic, with 
algebraic order of 0i. As the temperature increases, the 
stiffness J defined as J/2 = j^^/T decreases from infin- 
ity according to microscopic calculation. When J reaches 
2/tt, the system undergoes a Kosterlitz-Thouless transi- 
tion and the algebraic order of 0i is destroyed by prolifer- 
ation of stripe dislocations, similar to the destruction of 
superfluidity by vortices in the XY model. The system 
becomes nematic with a short-range correlation of 0i, 
while the C4 rotational symmetry remains broken. Fur- 
ther increasing of temperature eventually drives a second- 
order Ising-ncmatic phase-transition, above which the C4 
rotational symmetry is restored with (|0ip) = (|02p), 
and the system becomes normal. 

For the commensurate case with momentum 2fci? = 
27rp'/p, where p' and p are relatively prime integers, an 
additional term wcos{p(f) is allowed in Equation ([TT|). 
The U(l) symmetry of 0i is reduced to Zp here, which 
means the action is invariant under translation ip — > 
(p + 2'Kp" /p, with non-negative integer p" < p. With 
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FIG. 5. (Color online) Phase diagram at 1/5 filling where 
p = 5. Red (dark gray) dashed line: determining the 
smectic— algebraic smectic phase-transition temperature at 
J — p^/Stt. Green (light gray) dashed line: determining the 
algebraic smectic— nematic phase-transition temperature at 
J = 2/tv. 



this additional cosine term, the theory of (p natu- 
rally reduces to the Zp compact clock model. Ac- 
cording to the renormalization-group analysis of the 
compact clock model [1^, our system undergoes the 
smectic— nematic— normal transition as one increases the 
temperature when 1 < p < 4, and the smectic— algebraic 
smectic— nematic— normal transition when p > 4. As 
p ^ 00 where the system approaches incommensurate, 
the smectic— algebraic smectic phase-transition tempera- 
ture reduces to zero, which is consistent with the incom- 
mensurate case discussed before. 

As a specific, nontrivial example to make connections 
between the field theory given by Equation pT|) and the 
microscopic model given by Equation ([1} , a system with 
commensurate 1/5 filling {p = 5) is studied. We focus 
on the temperature regime near the Ising-nematic phase- 
transition point, so that in Equation (|lll) the approxima- 
tion — — r/2u for the saddle point is applicable. The 
procedure is to write Equation ([Ij in path-integral form 
and use Hubbard-Stratonovich transformation to intro- 
duce the density fields and decouple the interaction term. 
The fermionic fields are then integrated out and the re- 
maining density fields reproduce the field theory given by 
Equation (fTT|) with coefficients calculated. The details of 
the procedure are shown in Appendix |B] After obtaining 
the coefficients r,j, u, and the phase-transition temper- 
ature can be determined from the value of J according 
to the previous discussion. The phase diagram is shown 
in Figure [5l 

Our system has the following advantages to study 
liquid-crystal phases. (1) There are no other compet- 
ing orders. A similar spinful condensed-matter system 
was studied [2^ , which showed that a stripe CDW order 
of momentum {2kF,2kF) is competing with a checker- 
board CDW order of coexisting momenta (2fci?,7r) and 
{TT,2kp) [s^l- Also, a similar spinful ultracold-atomic 
system that has much more complicated interaction [iTj 
may have many competing phases such as BCS, which is 
not as clean and simple as the proposed spinless system 
to study the liquid-crystal phases. (2) The 2kp momen- 
tum dependence of density wave order is highly tunable 
by changing the fillings, which makes it easy to adjust 
the commensurability. 
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VI. EXPERIMENTAL REALIZATION 

Our system can be realized by loading fermionic atoms 
such as '"^K or ®Li of a single hyperfine state on square 
optical lattices. The fermions will automatically occupy 
p-orbital bands after the s band is fully filled. The in- 
teractions between spinless fermions can be tuned by p- 
wave Feshbach resonance, together with controlling the 
lattice spacing and potential depth. The cold gas is not 
required to be so close to the resonance, so the atom loss 
rate can be kept relatively low. The momenta of den- 
sity wave orders Qi,2 can be detected by optical Bragg 
scattering |3l| . Alternatively, in situ imaging can di- 
rectly show the density pattern Sp. As discussed be- 
fore, at half filling the checkerboard density waves are 
greatly enhanced by the Umklapp process, and are not 
affected by the Fermi-surface curvature caused by trans- 
verse hopping. Therefore, in experiments one should first 
search for the checkerboard density waves at half filling, 
which has the phase-transition temperature Tc ~ 0.2i|| 
based on RPA (a mean-field-level estimate) as shown in 
Sec. Illli given that = 2. In addition, the hopping 

t|| of p-band fermions is in general an order of magnitude 
larger than the s-band hopping, which also enhances the 
phase-transition temperature. With typical experimen- 
tal parameters such that A ^ 500nm, Vx ~ 5i?fl,, and 
using ^''K atom, the estimated parallel hopping from 
harmonic approximation is t|| ^ lOOnk, and then the 
phase-transition temperatures in Figure |3] and [5] can be 
determined. For example, the estimated phase-transition 
temperature for the checkerboard density waves at half 
filling is Tc - 20 nk from Figured 
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Appendix A: Validity of RPA 

We justify the validity of RPA used in Sec. [ITT] as fol- 
lows. Our system can be viewed as two sets of non- 
interacting one-dimensional (ID) spinless Fermi chains 
perpendicular to each other, with one set in the x- 
direction and the other in the y-direction. The x chains 
{y chains) are weakly coupled by small transverse inter- 
chain hopping in the y (x) direction, and show quasi-one- 
dimensional Fermi surfaces (Figure [1]). The interorbital 
interaction between and Py fermions is then turned on, 
which couples the motion of particles in the x direction 
and that in the y direction. In general, weakly coupled 
ID Fermi chains with intrachain interaction cannot be 
studied by RPA because of the Luttinger liquid behaviors 



in such quasi-one-dimcnsional systems. The well-defined 
(fermionic) single-particle excitations, which are required 
by RPA, are absent in Luttinger liquids, which makes 
RPA invalid in such weakly coupled ID Fermi chains [s^ . 
However, in our system, the key difference is the exis- 
tence of the interorbital interaction that couples p^ and 
Py fermions, which makes our system two dimensional 
(2D). This can be understood as follows. (1) An elastic- 
scattering process in ID between two particles with equal 
mass cannot change the momentum distribution of the 
two particles. Therefore, in a ID chain or weakly coupled 
ID chains with intrachain interaction, the momenta dis- 
tribution A^(k) = (C'''(k)C(k)) cannot be changed by the 
elastic-scattering process, which causes the excitations to 
be collective and the system to be a Luttinger liquid, for 
which the RPA is invalid. In contrast, in our system, the 
interorbital scattering process is a 2D scattering process, 
which can change the momentum and energy of the p^ 
fcrmion in x chains by transferring some momentum and 
energy to the Py fermion in y chains during the elastic- 
scattering process, and vice versa. In other words, the 
momentum distribution NxCk) = (C2(k)C2;(k)) of the 
Px fermions can be changed by scattering, and the same 
for the Py fermions. As a result, although our system ap- 
pears to be composed of ID chains of weak interchain 
tunneling, its dynamics is fundamentally 2D. At high 
temperature, a 2D system is at the Fermi-liquid phase 
with well-defined (fermionic) single-particle excitations, 
which is essentially different from the Luttinger liquid 
phase in ID or quasi-one-dimensional systems. There- 
fore, RPA is valid in our 2D system. (2) We can also 
classify the interactions according to their relevance in 
the renormalization-group (RG) flow to the Fermi sur- 
faces psl [ssj . The dominant marginal term in our system 
is 5Ek„k. CUk^ + Q)C.(ki)Ct(k2 - Q)C,(k2), where 
ki, ko, ki -I- Q and k2 — Q are all near the Fermi sur- 
faces [25| . This dominant interspecies interaction also in- 
duces effective intraspccics interactions. The induced in- 
traspecics forward and backward scattering processes are 
also marginal, but they are much weaker since they are 
higher-order processes. With the effect of the intraspecies 
interaction on the interspecies interaction neglected, the 
RG flow of the interspecies interaction diverges for tem- 
perature T < < exp(— ait/l^l) {ai is some constant). We 
emphasize here that this interaction describes scattering 
processes in 2D rather than ID. For such a 2D system 
with a single dominant interaction term, the RPA is well 
justified [2a |. 



Appendix B: Derivation of the Field Theory 

The coupling constants j, r, u, and v in Equation pip 
of Sec. |V] can be derived from the microscopic model. 
The Fcrmi-Hubbard model given by Equation ([T]) can be 
written in path integral form, where the partition func- 
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tion is given by e and the effective action reads 
Sf= j dr^?/>*(r,T)(9^ - /^)'0„(r,T) 

+ X!*"/3W'a(i" + e;3,T)V'a(r, t) + h.c.) 

ro/3 

+5 51 T)V'y (r, T)V'.(r, r). (Bl) 

r 

In Equation (|Bip . the interaction term can be rewritten 
as 

(B2) 

Consider two auxiliary Hubbard-Stratonovich fields 
/i?(pi,2)e^- ^ where = / drf P?,2(r, r). By 

shifting pi_2 -> Pi, 2 - (^^^x ± Pi and p2 denote 

the total density field and density difference field, respec- 
tively. Multiplying / D{pi^2)e^''^'''- with shifted pi, 2 to 
e~'^^ , the quartic interaction between fermions in Equa- 
tion (|B2p is eliminated. According to the mean- field 
analysis of CDW, the density difference has the saddle 
point (thermal average) at zero, which means we can ig- 
nore the density difference field p2, since the fluctuation 
around zero is trivial. As a result, the interaction term 
in Eq. ()B1|) is replaced by 



(B3) 



From mean-field analysis, the total density pi is fluctuat- 
ing around momenta 0, ±Qi, and ±Q2. The fluctuation 
around zero momentum is background density fluctua- 
tion, which is trivial. By ignoring such contribution, pi 
reduces to 5p, which reproduces the total density fluctu- 
ation 5p around Qi,2 in Sec. |Vl In the long- wavelength 
limit, the density fluctuations around Qi,2 can be rewrit- 
ten as 

+C.C.], (B4) 

by Fourier transform, where A is some momentum cutoff 
of the long- wavelength limit, and a — 1,2. Recall that 



in Sec. |Vl the flelds are defined through 

<5p(r,r)=^[0,(r,r)e^Q-'- + c.c.]. (B5) 



By Fourier transform. 



0<t(i-,t) 



T 



|<A,i 



=i(q-i 



^)0,(q,w), (B6) 



and comparing Equation (jB4p with Equation (jB5| . we 
reach the relationship (5/c(Qo- + q, = ^^(q, w). The ef- 
fective action is then written in momentum space, where 
5 Pa can be replaced by ^g-. At last, the fermionic fields 
are integrated out and the (j) fields are kept up to quartic 
terms, and we reach the expression 



+t;|0ip|02|2 + 5t + .... 



(B7) 



which reproduces the field theory given by Equation (jlip 
in Sec. |Vl The coefficients in Equation (|B7|) are 
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Jl 

' 2T m^T 



E 



1 ~ 2nf (Ck) 
2Ck ' 



16iV2T 



E 



d'^np sinfc^.)^ 
1 — 2nF dnp 1 



32iV2T ^ 

k 



V = 4u, 



(B8) 



in static limit (we do not consider quantum fiuctuations 
in the present work). Here, ^k is the spectrum of free 
Px-orbital fermions. 

In Sec. |Vl a 1/5 commensurate filling case is consid- 
ered. In order to obtain the coefficients in Eq. (jB8|) in this 
case by the above procedure, the chemical potential is ad- 
justed to make the filling 1/5 {kp = 4/57r), with CDW 
momentum Q = (^, ^). A term wcos(5</j) produced by 
+ C.C. also exists, w in general is small because this 
term arises from high-order diagrams and is suppressed 
at flnite temperature. The exact value of w is not im- 
portant, but the existence of this term is crucial in the 
commensurate filling. The field theory then reproduces 
the Zp compact clock model, and the RG analysis of this 
model can be used to determine the phase-transition tem- 
perature, as discussed in Sec. [V] 
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